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We present a calculation of the electromagnetic form factor of the pion in Nf = 2 + 1 flavor 
lattice QCD. Calculations are made on the PACS-CS gauge field configurations generated 
using Iwasaki gauge action and Wilson-clover quark action on a 32"^ x 64 lattice volume with 
the lattice spacing estimated as a = 0.0907(13) fm at the physical point. Measurements of the 
form factor are made using the technique of partially twisted boundary condition to reach 
small momentum transfer as well as periodic boundary condition with integer momenta. 
Additional improvements including random wall source techniques and a judicious choice 
of momenta carried by the incoming and outgoing quarks are employed for error reduction. 
Analyzing the form factor data for the pion mass at ~ 411 MeV and 296 MeV, we 
find that the NNLO SU(2) chiral perturbation theory fit yields (r^) = 0.441 ± 0.046 fm^ 
for the pion charge radius at the physical pion mass. Albeit the error is quite large, this is 
consistent with the experimental value of 0.452 ± 0.011 fm^. Below w 300 MeV, we find 
that statistical fluctuations in the pion two- and three-point functions become too large to 
extract statistically meaningful averages on a 32"^ spatial volume. We carry out a sample 
calculation on a 64'' lattice with the quark masses close to the physical point, which suggests 
that form factor calculations at the physical point become feasible by enlarging lattice sizes 
to M^L « 4. 



I. INTRODUCTION 



The electromagnetic form factor of pion is an interesting quantity to investigate in lattice 
QCD. Experimentally it has been measured in a set of experiments [l[ . Together with the nucleon 
form factors, it provides the first important test case of our understanding of hadron structure 
that arises from the quark content. From lattice QCD point of view, form factor calculations 
represent one of the first steps going beyond static quantities like the mass spectrum which 
require only two-point functions. The pion form factor is a natural first choice in this direction 
since usually pion Green's functions are statistically the most stable quantities in lattice QCD 
measurements. An interesting point with the pion form factor G-j^^q^) is its slope at the origin 
as a function of the momentum transfer squared q^, i.e., the pion charge radius defined by 

It has been known for some time from chiral perturbation theory analysis that this quantity 
diverges logarithmically with vanishing pion mass squared. Quantitative confirmation of such 
a behavior would provide an important check on the control of chiral behavior in lattice QCD 
simulations toward the physical point. 

The pioneering lattice QCD calculations of the pion form factor appeared more than 20 
years ago 0, Q], and a number of studies were carried out over the years. Recently, with the 
development of simulations with dynamical quarks, several groups have attempted calculations 
with A'^^ = 2 0-0] and Nj = 2 -|- 1 [o^ dynamical flavors using a variety of quark actions. The 
Nf = 2 calculations employed Wilson-clover [6], twisted mass {3] or overlap jsl quark action, 
and explored the pion mass region down to M^^ ~ 300 MeV. The values for (r-^) from those 
simulations are significantly smaller than the experimental value, and NNLO fits of SU(2) chiral 
perturbation theory were needed to find consistency with it at the physical pion mass. For 
Nf = 2 -|- 1 dynamical fiavors, there has been a single calculation employing domain-wall quark 
action [9], which made measurements at a single pion mass of ~ 300 MeV. Carrying out NLO 
analyses in SU(2) and SU(3) chiral perturbation theory, this work found (r^) to be consistent 
with the experimental value at the upper edge of a 10% error band. 

In this paper we present our calculation of the electromagnetic form factor of pion in Nf = 
2 + 1 dynamical flavor QCD using the Wilson-clover quark action. For measurements we employ 
the Nf = 2 -|- 1 PACS-CS gauge configurations generated on a 32^ X 64 lattice using the Iwasaki 
gauge action and the Wilson-clover action at a lattice spacing estimated to be a = 0.0907(13) fm 



at the physical point [10]. Since the pion mass on this gauge configuration set covers the range 
from ~ 700MeV down to 156 MeV, we are able to examine both the known range above 
M,r ~ 300MeV and a novel range below toward the physical pion mass. 

The paper is organized as follows. In Sec. 2 we present our method to calculate the pion 
form factor. In order to access the region of small momentum transfer, we use the method 



of partially twisted boundary condition 



11 



14l |. and in order to fight increasing computational 



a, 



15l4l9l|. In addition 



cost for smaller pion mass, we apply the method of random wall source 
we make use of a judicious choice of momenta carried by the incoming and outgoing quarks off 
the electromagnetic vertex, which helps in reducing statistical fluctuations in the form factor 
measurements. In Sec. 3 we present the results of pion form factor measurements, and in Sec. 4 
we analyze the data as a function of the momentum transfer squared and pion mass. In particular 
we examine consistency with the predictions of chiral perturbation theory. Finally, in Sec. 5, we 
discuss our findings closer to the physical point including the results of our test calculation on 
a 64^ lattice with the quark masses tuned to the neighbour of the physical point. We end this 
work with conclusions in Sec. 6. A preliminary report of this work was presented in {2^. 



II. METHODS 



A. Pion electromagnetic form factor 



The electromagnetic pion form factor Gt^{Q'^) is defined by 

J^|7r+(p-)) = {p^+p'^)GAQ^), (2) 

where = —q^ = —{p' — p)"^ is the four-momentum transfer, and is the electromagnetic 
current given in Nj = 2 + 1 QCD by 

2_ 1 - 1_ 

Jti = -^u-i^u - -d'y^d - -S7^s. (3) 

In the limit of vanishing four-momentum transfer = 0, the form factor equals unity due to 
the charge conservation. 

We extract the form factor from a suitable ratio of the pion two- and three-point functions. 
We use a ratio, which has the advantage of simultaneously reducing fluctuations and renormal- 
izing the current, defined as 

^ C^P\p',tf;p,0;T) C^p\p',t) ^ 2E^{p>) 
C^P^{p',tf;p',0;T)C^P'ip,r) E^{p) + E^{p>y 



which converges as 

^ = GAQ% (5) 

for large r and tj. Et^{p) denotes pion energy for spatial momentum p, C'^^^{p' ,tf;p, 0; r) is the 
three-point function with momenta p at the source and p' at the sink tf, 

C^P\p',tf;p,0;T) = {7r+{p',tf)J^{T)7T+{p,0)) , (6) 

and C'^p^{p,t) is the two-point function, 

C2P*(p,r) = (7r+(p,r)7r+(p,0)). (7) 

After contraction of quark fields, the three-point function consists of the connected and discon- 
nected contributions. The latter vanishes after the gauge field average due to charge conjugation 
invariance, and hence need not be calculated. Since we assume degeneracy of up and down quark 
masses in the present calculation, the connected contribution is equal to 



C^P\p',tf,pAr) =Y.e-^P'y+^'^^{Tr[roD-\Q,x)^,D-\x,y)^,D-\yM) ■ (8) 



This contribution can be calculated by the traditional source method 



B. Choice of momenta carried by quarks 



The ratio dH) makes use of two- and three-point functions in an appropriate combination 
to extract the form factor for the renormalized current. The presence of ratios guarantee that 
statistical fluctuations are suppressed. Nonetheless, making simple choices such as p ^ and 
p' = 0, we have observed an increasingly larger fluctuation of the ratio as pion mass is reduced, 
and this trend worsens for larger momenta. With an interesting choice of momenta, p' ^ p but 
\p'\ = \p\, the ratio @ simplifies to 

C!»"(p',t/;p',0;T) 

Since the two-point functions as well as the ratio of energies drop out, leaving just the ratio of 
three-point functions, we expect this choice to yield better signals than those choices for which 
all factors are present. Furthermore, one can choose 6 permutations in momentum directions 
while keeping \p'\ = \p\, gaining more statistics. 
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C. Partially twisted boundary condition 



The minimum non-zero quark momentum 2tt/ La for the periodic boundary condition on a 
32^ X 64 lattice with a 2 GeV inverse lattice spacing is about 0.4 GeV. To probe the region of 
smaller momentum transfer as well as to improve the resolution of four- momentum transfer, we 
apply the method of partially twisted boundary condition Ill4l4i| in which valence quark fields 



are subjected to twisted boundary condition while periodic boundary condition is kept for sea 
quark fields. If one imposes the boundary condition given by 

ip{x + Lej) = e^''''^^i^{x), j = 1,2,2,, (10) 

on a valence quark field, the spatial momentum of that quark is quantized according to 

Pi = -^ + -^, j = 1,2,3, (11) 

where L denotes the spatial lattice size, Cj the unit vector in the spatial j-th direction and 9j real 
parameter. In this way one can explore arbitrarily small momentum on the lattice by adjusting 
the value of twist 9j. 

For the meson two-point function consisting of quark and anti-quark propagators, we apply 
the twist only to quark and not to antiquark or vice versa. Similarly, for the three-point function, 
we twist only one or two out of the three quark propagators. In other words, we pretend that 
each valence quark line in the two- and three-point function quark diagrams carry a different 
flavor and select the appropriate flavor to apply twisting. In this way we can avoid a twist 
of a quark line cancelled by the opposite twist of the antiquark line carrying the same flavor 



lll |. This procedure and the twisting of only valence quarks mean that we deal with partially 



quenched QCD with a different flavor symmetry content in the valence and sea quark sectors. 
As was discussed in detail in [l^, 3] using chiral perturbation theory, the associated effects are 
expected to appear as finite-size effects exponentially small in spatial volume for channels which 
do not have final-state interactions such as three-point functions for form factor calculations. 
Since terms of such magnitude are also present in unitary theory with periodic valence and sea 
quarks, we ignore this issue in the present work. 

The twisted boundary condition can be imposed on a periodic quark field configuration by 
the following transformation 

il){x) U{e,x)il){x) = e^"^i=i^^"'^/'^V(a;). (12) 




_j 1 1 1 1 1 — I o_ii 1 1 1 1 1 1 1 1 1 

5 10 15 20 25 30 0-1 0-2 0-3 0.4 0-5 0.6 0.7 0. 



time 0^ 
(a)Pion effective energies. (b)Energy momentum dispersion relation. 

FIG. 1. Check of validity of the twisted boundary condition at Ks = 0.1364, = 0.13700 where 
~ 702 MeV. Values of twist angle used in test are 9 — 0.1842 and 0.2811. Measurements are made 
on 40 configurations. 



In practice we transfer the twist from the quark sector to the gluon sector by an U(l) transfor- 
mation on the spatial gluon Hnk fields given by 

Uiix)^u!{x) = e^'''^^/^Uiix), i = 1,2,3. (13) 

Thus valence quark propagators are solved with the periodic boundary condition but on the 
PACS-CS gauge configurations twisted by the U(l) transformation above. 

In order to check that the term acts as true physical momentum, we carried out a test of 
the energy-momentum dispersion relation of the pion on some PACS-CS configurations. Of the 
two valence quarks inside the pion, we twisted one quark with a twist angle 9 = {9, 9, 9) and left 
the other untwisted. In Fig. HJa) we plot the effective energy for the ground state, two values of 
the twist angle and their combination with the first integer momenta at the hopping parameters 

= 0.1364, Kud = 0.13700 where 702 MeV. The propagator is fitted over t = 7 - 27 to 

extract the energy E{p). Errors are estimated by the jackknife method with the bin size of 100 
trajectories. The results are plotted in Fig. [Hb), together with the expected behavior, 

E{p)' = iaM^r+(^^n+^-^9y, (14) 

which demonstrates clearly that the term acts as true physical momentum. The two data 
points on the right represent combinations of an integer momentum (1,0,0) and a twist. The 
energy-momentum relation is correctly reproduced in this case as well. 
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D. Random wall source 



At light quark masses the computing cost for inversion of Dirac operator becomes very 
expensive. Thus we have employed some improvements for obtaining the form factor with 
acceptable statistical errors at reasonable computing time. The first improvement is to utilize 
the random wall source. This method has a long history and has been applied to two-point 
functions in a variety of contexts. More recently, applications to three-point functions have 
shown their effectiveness for form factor calculations [3, We consider the use of Z{2) Z{2) 
random noisy source as introduced in [17]. 

Consider a set of random sources whose real and imaginary components are randomly chosen 
from Z{2) for each site, color and spin, 

{r/(")(3;)aa G Z{2) ® Z{2)\n = I...N}. (15) 

This set has the property that in the limit N ^ oo 

1 N 

n=l 

To use this kind of source in calculating correlators, one can choose it to be a set of random 
wall source located at to, 

r?(^)(x,t|to) eZ(2)$$Z(2)|t = to 

= 0\t + to, (17) 

{vlu^{x,t\to)vl''^\y,t\to))n = 6^y6ab6a/3, N ^ oo. (18) 
Making use of (jlSp to rewrite the pion two-point function at zero momentum as, 



C(t; 0) = E tr [D-^y, t; x, to)D-^\y, t; x, to) 

x,y 

= {^aaMy^ *o) [SxzSbc^tsK] D^K,L{y, z, to) 

x,y,z 

= E(^^"^(y'*IW^^"^(y,i|to))^, (19) 

y 

where ip^"'^ is the solution vector of the Dirac equation, 

i^^''\y,t\to) = J2D-\y,t;x,to)r]^''\x,t\to). (20) 
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(a) Effective pion mass with point source. 
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(b)Effective pion mass with random Z{2) (g) Z{2) wall 
source. 
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(c)Effective pion mass with smeared source. 

FIG. 2. Comparison of pion effective masses calculated with (a) point source, (b) random Z{2) (E) Z{2) 
wall source, and (c) smeared source on a set of 10 configurations at Kg ~ 0.1364, K^^i = 0.13700 (M^r « 
702 MeV). 



With a random Z(2) <^ Z{2) wall source the solution for quark propagator needs only single 
inversion instead of 3 x 4 = 12 corresponding to color and Dirac components required for a point 
source or smeared source. When the number of configurations in the ensemble is large enough, 
even if one uses a single random source for each configuration, (llSp is expected to hold in the 
ensemble average. One may expect to obtain meson correlators of a similar statistical quality 
as with the traditional point source with only 1/12 of computing time. 

In Fig. [21 we compare the effective pion mass plot calculated on a set of 10 configurations 
at Ks = 0.1364, = 0.13700 (M^ ^ 702 MeV) from the PACS-CS ensemble using (a) point 
source, (b) single random Z(2) Z{2) wall source, and (c) smeared source. We observe that 
the signal with the random Z{2) ® Z{2) wall source is somewhat better than that for point 
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source, while the signal for smeared source is better than that with the random Z{2) (8) Z{2) 
wall source. Using 4 random wall sources for each configuration, we observed that the quality 
of signal becomes comparable to that for smeared source. Since the computing time is still 
4/12 = 1/3 for the random wall source, we employ the method of random wall source with 4 
sets of random wall sources in our measurements. In addition we repeat measurements with the 
source located at i = 0, 16, 32, 48 since the time extent of our lattice is 64. 



We apply our calculational setup to a subset of the PACS-CS gauge configurations [1; 



corresponding to the degenerate up-down hopping parameter in the set Kud — {0.13700, 0.13727, 
0.13754, 0.13770}. The hopping parameter of strange quark is fixed at Kg = 0.1364. 

The first set of measurements, which we call data set I, is made with an exponentially smeared 
source and local sink, setting the final pion at zero momentum p' = and varying that of the 
initial pion p in the three-point function. The fixed sink time tf in the ratio ^ needs to be 
chosen large enough to eliminate excited states contributions. However, statistical fluctuations 
increase as tf increases, and examining measurement results, we choose = 24 to balance the 
two opposite features. The twist technique is applied to the quark running from the source 
to the current. Two values are chosen for the twist angle 6 = {9,6,9) such that the smallest 
four-momentum transfer of the current takes the value Q^(GeV^) = 0.01841 or 0.04237. Adding 
integer momenta, we then collect data for in the range 0.01841 GeV^ < < 0.7302 GeV^. 
The statistics of data set I is given in Table H] together with pion and kaon mass. Results of data 
set I have been previously reported in {2^. 

In order to extract the form factor, we fit the plateau of the ratio R{t) by a constant. The 
fitting range should be chosen around the symmetry point between the source and the sink, with 
additional considerations on the time interval required for the pion state to become dominant. 
Since we employ an exponential smeared source and a point sink, we shift the fitting range one 
time unit closer to the source than the symmetric point tf/2 = 12. 

In Fig. [3|^a) we plot the ratio R{t) at various momentum transfer for the pion mass ~ 
702 MeV. At this pion mass we have good signals for all 7 values of the four-momentum transfer. 
There is a good plateau from r = 8 to 15 for every momentum transfer. Thus at this pion mass 
we can choose the fitting range from r = 8 to 15 to extract the form factor. However, as 



III. MEASUREMENTS 
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(a)Data set I: Ratio R{t) defined in Q at all 7 values (b)Data set II: Ratio R'{t) defined in ^ at all 5 values 
of 4-momentum transfer. of 4-momentum transfer. 

FIG. 3. Ratios to extract the form factor as functions of the time shce t of the current operator at 
A'U ~ 702 MeV. 



(MeV) Mk (MeV) #conf measured 



0.13700 0.1364 


702 


792 


40 


0.18423, 0.28112 


0.13727 0.1364 


570 


716 


40 


0.18467, 0.28265 


0.13754 0.1364 


411 


637 


40 


0.18585, 0.28672 


0.13770 0.1364 


296 


596 


160 


0.18814, 0.29450 



TABLE I. Statistics of data set I. 

the pion mass decreases, the plateau signal becomes worse as exhibited in Fig. ID^a) for the 
lightest case of M,r ~ 296 MeV where the ratio R{t) at the two smallest momentum transfers, 
(5^(GeV^) = 0.01841 and 0.04237, is shown. We then choose larger values for the starting 
point of the fitting range for better suppression of excited states at lighter pion masses. The 
error is estimated by the jackknife method using 10 configurations corresponding to 50 hybrid 
molecular dynamics time units as the bin size after checking saturation of the magnitude of error 
as function of bin size. Fit results for the pion form factor are listed in Table HIl 

We observe in Table HH for data set I that the error for the form factor becomes large toward 
small pion mass and large momentum transfer. In order to improve the quality of data, we 
repeat the measurements (i) choosing the incoming and outgoing pions to have momenta with 
the same magnitude \p'\ = \p\, and (ii) applying 4 random Z{2) ® Z{2) wall sources located 
at t = 0, 16, 32, 48 for the lattice time extent of 64. The twist technique is applied to two 
quarks running from the source to the current and from the current to the sink. Five values are 
employed for the twist angle of form 9 = (0, 0, 0) and its permutations such that four- momentum 
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(a)Data set I: Ratio R{t) defined in Q at 2 smallest (b)Data set II: Ratio R'{t) defined in ((9]) at 4 smallest 
4-momentum transfers: Q^{GeV^) = 0.01841,0.04237. 4-momentum transfers: 

g^(GeV^) = 0.02, 0.04, 0.06, 0.08. 

FIG. 4. Ratios to extract the form factor as functions of the time slice t of the current operator at 
AU ~ 296 MeV. 



Mrr ~ 702 MeV, fit range: [8-15], bin size: 50t 
Q2(GeV2) 0.01842 0.04237 0.1163 0.1258 0.1682 0.3651 0.7302 
G^(Q2) .9825(24) .9609(43) .8834(120) .8780(134) .8511(188) .7313(186) .5875(200) 



AU ~ 570 MeV, fit range: [9-14], bin size: 50t 
Q2(GeV2) 0.01842 0.04237 0.1132 0.1223 0.1623 0.3651 0.7302 
G^(Q2) .9836(37) .9604(61) .8816(154) .8746(160) .8400(184) .6934(212) .5191(189) 



~ 411 MeV, fit range:[9-14], bin size: 50t 
Q2(GeV2) 0.01841 0.04237 0.1062 0.1143 0.1495 0.3651 0.7302 
G^(Q2) .9730(54) .9428(66) .9036(315) .8920(319) .8805(476) .5999(535) .4706(574) 



« 296 MeV, fit range: [10-13], bin size: 50r, 
Q2(GeV2) 0.01842 0.04237 0.09612 0.1030 0.1324 0.3651 0.7302 
G^(Q2) .9728(44) .9372(72) .8624(310) .8456(343) .7929(452) .9758(3376) .6115(1943) 

TABLE II. Pion form factor obtained with data set I. 

transfer of the current takes the value Q'^{GeY^) = 0.02, 0.04, 0.06, 0.08, 0.10. Those values of 6 
are independent of as is easily checked for the momentum configuration chosen here. The 
fixed sink time is chosen to be 28, larger than that of data set I, for better suppression of 
excited states and also from examination of the dependence of the ratio on tf. The fitting range 
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is chosen symmetric around t = 14 since the source is local after averaging over the Z(2) (g) Z{2) 
random numbers. We call this set of data as data set II. Statistics and results of data set II 
are tabulated in Table Hill and Table ITVl Results for R'{t) for M-j^ w 702 MeV are plotted in 
Fig. mb). One can see that the form factors of the data set II have much smaller error bars 
compared to those of data set I. We also plot results for the case of pion mass ^ 296MeV 
in Fig. ifb). 

Kud Us Mtt (MeV) Mk (MeV) #conf measured GeV 2) 

0.13700 0.1364 702 792 40 0.02, 0.04, 0.06, 0.08, 1.0 

0.13727 0.1364 570 716 40 0.02, 0.04, 0.06, 0.08, 1.0 

0.13754 0.1364 411 637 40 0.02, 0.04, 0.06, 0.08, 1.0 

0.13770 0.1364 296 596 160 0.02, 0.04, 0.06, 0.08, 1.0 

TABLE III. Statistics of data set II. 



Q2 (GeV2), fit range: [11,17] 



M^(MeV) bin size 0.02 0.04 0.06 0.08 0.10 

702 50t 0.9818(5) 0.9645(9) 0.9473(17) 0.9308(17) 0.9155(22) 
570 50r 0.9796(6) 0.9562(17) 0.9385(23) 0.9217(27) 0.9030(31) 
411 50t 0.9727(11) 0.9506(23) 0.9229(34) 0.9083(52) 0.8927(75) 
296 50t 0.9733(16) 0.9462(45) 0.9221(50) 0.8911(70) 0.8959(96) 

TABLE IV. Pion form factor obtained with data set II 



IV. PION ELECTROMAGNETIC FORM FACTOR AND CHARGE RADIUS 
A. Monopole analysis of the dependence of the form factor 

Figure [5] shows the momentum transfer dependence of our data for the pion form factor 
at all simulated pion masses. The data set I and II are consistent with each other within the 
estimated errors. The experimental pion form factor is phenomenologically reasonably described 
by a monopole form suggested by the vector meson dominance model, 



Our data are accordant with the ansatz; solid lines in Fig. [5] are fits to the monopole form ()2ip . 
For monopole analysis, we utilize the form factor data in the range up to = 0.08 GeV^ at 
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FIG. 5. Momentum dependence of Gt^{Q^) from data set I and II. Blue lines are fits of data I and II to 
the monopole ansatz ([21]). 



= 296 MeV and up to 0.10 GeV'^ at 411 MeV, since at larger four-momentum transfers 
plateau signals are not clear. The fitted values of the monopole mass M^^^^ can be used to 
estimate the pion electromagnetic charge radius via (r^) = G/M^^^^. Results are tabulated in 
Table IVl and plotted in Fig. [6l The charge radius exhibits an increase as pion mass decreases. 



Kud 0.13700 0.13727 0.13754 0.13770 
(r2)(fm2) 0.2174(27) 0.2538(38) 0.3129(84) 0.3352(160) 
X^/d.o.f 0.51(12) 1.01(12) 1.12(19) 0.44(11) 

TABLE V. Pion squared charge radius calculated from the monopole fit of data set I and II. Errors are 
estimated by Jackknife method with bin size of 50 t. 
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► This work 

* Experimental value 



M^(GeV^) 



FIG. 6. Squared charge radius (r^^ [br?) as a function of A/^ obtained with the monopole ansatz. The 
left-most point represents experiment. 

B. ChPT analysis to NLO 

For small values of momentum transfer and pseudo-scalar meson masses, we expect ChPT to 
provide a description of the pion form factor as a function of those variables. Here we analyze 
our data in terms of ChPT to NLO. The analytical expression for the form factor has been 
worked out long time ago both for SU(2)l x SU(2)r and SU(3)l x SU(3)r [s] cases, which 
is given by 



2Ml r sL H{s) 
T ^ N 



and 



4Q^ 



2Mi 



sL His) 

\ — 

6 N 



+ 



Ml 



srLk . H{sk) 



6 



N 



(22) 



(23) 



where 



H{x) 



x~4 



4 5 X - 4 /x - 4, 

1 X \ log 

3 18 6 V X ^ \ Ix-i 



+ 1 



(24) 



and / and /q are the decay constant in the SU(2) and SU(3) chiral limit, respectively. In the 
above equations, we made use of the following definitions: 



= (47r)^ 
Ml 



(25) 



■■,SK 
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M^(MeV) SU(2) SU(3) 

- l/«) x'/d.o.f {r')^^Jfm^) = 1/a) xVd.o.f {r')^^Jfm^) 

296 -0.00737(45) 0.29(10) 0.366(14) 0.00256(19) 0.29(10) 0.380(14) 

411 -0.00728(26) 2.34(26) 0.363(8) 0.00260(11) 2.31(26) 0.383(8) 

TABLE VI. NLO ChPT fit of G^{Q^) at fixed pion mass AU. LEC's are calculated at fi ^ 1/a = 
2.176 GeV. The decay constants in the chiral limit are taken from the work of PACS-CS collaboration 
Q: / = 124.8(5.1) MeV,/o = 116.0(8.8) MeV. Those values are determined with pion masses up to 
411 MeV. Physical value of the squared charge radius, {f'^)p}^ygi is calculated at the physical pion mass 
for the SU(2) case and at the physical pion and kaon masses for the SU(3) case. 

Besides the decay constant at the chiral hmit, the SU(2) formula (|22p involves 1^ as the only 
unknown LEG, and the same situation holds for the SU(3) case, (j23p . with Lg as the unknown 
LEG. Galculating the slope at the origin of the momentum transfer yields the expressions for 
the charge radius: 

In Table rvTl we present results of NLO fits of the form factor for both SU(2) and SU(3) GhPT. 
Fits are made fixing the pion mass at each of the two lightest values, M-,^ = 411 and 296 MeV, 
available in our data set. Let us first look at the SU(3) results. In this case the measured kaon 
and pion masses are used in the fit, while the physical masses are substituted for computing 
the charge radius at the physical point from the fit results. The charge radius extrapolated to 
the physical point, while consistent within the error for the two pion mass values, falls short of 
the experiment by about 15%. The SU(2) results in Table IVII are similar. The value for the 
charge radius predicted at the physical pion mass is about 20% smaller than experiment. We 
find similar values in the previous studies [g-^ carried out over a similar range of pion mass. 

In Fig. Owe plot the fit curves of the pion form factor together with curves from the monopole 
ansatz for the case of = 296 MeV for (a) SU(2) and (b) SU(3) GhPT to NLO. The pion mass 
dependence of the squared charge radius (r^) (fm^) which results from the fits are given in the 
panels (c) and (d) for the SU(2) and SU(3) cases, respectively. Filled circles are the estimates 
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(a)NLO SU(2) ChPT fit of the form factor. 



(b)NLO SU(3) ChPT fit of the form factor. 



* Experimental value 

• Monopole ansatz 
NLO SU(2) ChPT 

C, NLO SU(2) ChPT at the physical M 
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* Experimental value 




• Monopole ansatz 




i NLO SU(3) ChPT 




NLO SU(3) ChPT at the physical 
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(c)(r-2) determined from the NLO SU(2) ChPT fit of id){r^) determined from the NLO SU(3) ChPT fit of 
the form factor in comparison with values from the the form factor in comparison with values from the 
monopole ansatz and the experiment. monopole ansatz and the experiment. Along the ChPT 

line displayed in the figure, the kaon mass is fixed at its 
physical value. 

FIG. 7. Results of NLO analysis at the pion mass Af^ « 296MeV. Fit curves as compared to the form 
factor data and results for the charge radius are shown for the SU(2) case in panel (a) and (c), and for 
the SU(3) case in panel (b) and (d). 



from the monopole ansatz, and the asterisk on the left is the experimental value. As indicated 
in the panels (c) and (d), the NLO ChPT predictions for (r^) at 296 MeV are smaller than that 
obtained from the monopole ansatz. These differences are also visible in the panels (a) and (b) 
as indicated in the magnified region of small four-momentum transfers. 

The NLO ChPT fit at = 296 MeV has a smaller /d.o.f compared to that of monopole 
fit tabulated in Table |Vl This is not the case at the pion mass of 411 MeV, however, where 
/d.o.f of the NLO ChPT fit is significantly larger than that of the monopole fit. As we shall 
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M7r(MeV ) 


QJJfO\ 

oV[Z) 






= l/a) xVd.o.f {r^)^^Jim^) L^if, = 1/a) 


xVd.o.f {r')^,Jfm^) 


296 


-0.01238(66) 0.29(10) 0.457(18) 0.00577(33) 


0.29(10) 0.462(18) 


411 


-0.01408(45) 2.38(27) 0.502(12) 0.00666(22) 


2.37(27) 0.509(12) 



TABLE VII. NLO ChPT fit of Gj^iQ^) utilizing M^/f^ as the expansion parameter at fixed pion mass 
Mtt- Values of the decay constant at the simulation points are /tt = 151.7(2.7)MeV and 162.8(2.6) MeV 
at AU = 296MeV and 411 MeV, respectively. The physical decay constant, fP^y' = 132.7(5.5)MeV, is 
obtained from analysis of data with pion masses up to 411 MeV. Physical value of the squared charge 
radius, {f'^)pi^yg, is calculated at the physical pion mass for the SU(2) case and and at the physical pion 
and kaon masses for the SU(3) case. 

discuss below in Fig. [5|^a) for the SU(2) case, this is due to an upward curvature of the form 
factor data as increases to 0.08 and 0.10 GeV^. Higher order terms in need to be included 
in order to explain the behavior of our data for the form factor at 411 MeV. 

We now investigate the choice of the decay constant to be used in the NLO ChPT fit. 
With an uncertainty of order O(p^), the decay constant in (|22p and (|23|) can be chosen to be f-,^ 
measured at each pion mass. Table IVTIl shows results for NLO fits using M^/f^ as the expansion 
parameter. We observe a large difference in the results depending on whether one uses / or f.,^ 
for the SU(2) case and /o or /,r for the SU(3) case. The fit results for /g and L'g if one uses fir are 
larger than those of the fit using / and /o by 40 to 60%, which raises the values of {'f'^)pfiyg at the 
physical point by 20 to 30%. Predictions for the charge radius from these fits overestimate the 
experimental value while those employing / underestimate it. This uncertainty clearly indicates 
the importance of terms in the ChPT interpretation of our form factor data at the considered 
range of pion mass. 

Comparison of results using / and /jr for the SU(2) case at 411 MeV are made in Fig. [8j The 
left panels (a) and (c) are results obtained with / while those using /^r are shown in the right 
panels (b) and (d). In both cases, NLO ChPT fits do not explain the upward curvature of the 
form factor data at = q.OS and 0.10 GeV^. 

We should note that the SU(2) ChPT analysis requires tuning of the strange quark mass 
to the physical value, or alternatively, the dependence of the SU(2) LEC's on around its 
physical value has to be determined from data. For the Wilson-clover quark action, there is an 
additional subtlety that the strange quark mass, as defined via the PCAC relation, varies with 
changing up-down quark hopping parameter even if the strange quark hopping parameter is kept 
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* Experimental value 
Monopole ansatz 
NLO SU(2) ChPT using f 
^ NLO SU(2) ChPT at the physical M 



* Experimental value 

• Monopole ansatz 

o NLO SU(2) ChPT using f 
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0.1 0.12 0.14 0.16 0.18 



(GeV^) 



(GeV^) 



(c)^r^^ determined from the ChPT fit of the form (d)^r^^ determined from the ChPT fit of the form 

factor utiUzing /. factor utilzing 

FIG. 8. Results of SU(2) NLO analysis at the pion mass Mt^ w 411MeV. Fit curves as compared to the 
form factor data and results for the charge radius are shown for the case of using / p as the expansion 
parameter in panel (a) and (c), and for the case of using Al"^/ f^^ in panel (b) and (d). 



fixed. Our data taken for only one value of the strange quark hopping parameter, however, is 
not detailed enough to fully resolve the rUg dependence. We leave such a precise determination 
of the rus dependence for future work. 



C. ChPT analysis to NNLO 

The ChPT analysis to NLO presented in the previous subsection indicates that the NLO is 
not sufficient for the pion mass as large as M^^ ~ 300 — 400 MeV. Attempts have been made to 
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carry out fits to NNLO of CliPT 0,0, and we try this procedure here for the SU(2) case. 
The NNLO formula for the vector form factor is given by 22 1 



^5C/(2),iViVLO(g2)^;L + 2x2 



-{s-A)J{s) + s 



"'■fi r 



1 



6 18A^ 



+ 4x1 (Pf +4'0+O(^2), 



(28) 



where 



M. 



-f2 



X2 



P 

J TT 



(29) 



with the decay constant at the pion mass M,r, which is related to the decay constant / in 
the SU(2) chiral limit through 



fn = f 



.Ml 
P 



(30) 



at NLO of SU(2) ChPT. The two functions P^^ and U^^ are given by 



1111 

-ki + -k2 - —ki + -/C6 - II ( 211 + 



1 \ 23 L 5 



37 



^ — ki - —k2 + —kQ + — (l^ - -lo + -Ir 

12 24 24 9iV V 2 2 2 ^ 12 384 192A^ 



47 \ 

) + ""^^ 



(31) 



= J 



hli-s' + 4s) + If,{s^ - 4s) + hlis - 4) + hli-s' + 4.) 
3 3 o 

L(s^ + 8s - 48) H s + - 

36 ^ ' N \ W8 108 4 

+^K,(s) + ii^2(s) (^s^ - s + 4) + ^-Ks{s) (s-l)- Ik,{s), 



(32) 



and the integral functions J , Ki, K2, K3, K4 are defined by 



J 










z 


-4iV 













z 












K2 












8 




h 


K3 




Nzs' 


-1 


7r2(iVs) 


-1 ^2 




-(2iV2)-i 



(33) 



and 



1/1 1 1 - (7r2-6)s 

2^^ + 3^^+iV^ + 42iV^ 



(34) 
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exp. parameter M^(MeV) l^{p = 1/a) lO'^r^i 10V^2 X^/d-o.f (r^) 



phys 





296, 411 


-0.0098(11) 1.67(51) 1.04(60) 0.63(9) 0.420(31) 




296, 411 


-0.0103(18) 3.4(1.8) 3.3(1.6) 0.72(10) 0.441(44) 



TABLE VIII. NNLO ChPT SU(2) fit of G^(Q2) using data at 2 iightest pion masses. The result in tlie 
first row is obtained by using / — 124.8(5.1) MeV in the SU(2) chiral limit while that in the second 
row by substituting the measured values of f-^ . 

where 

h{s) = ^ln^^^, z = l--. (35) 
As well as notations in (j25p . we also use 

h = [4/[ - -fiL]L, (36) 

where 

71 = 1/3, 72 = 2/3, 74 = 2, 76 = -1/3. (37) 
From ()28p the NNLO ChPT expression for the squared charge radius reads, 

= -- 1 fe/^ + ^ + T?) (38) 

/ SU 2),NNLO fZ\^N 



SU(2),NNLO N 

4M2 r ^, 3, k4 ^, 1 / 31 13 181 \ 

— ^ -3ki + -k2 - + Ska - UrAl H -2V, H L H + QrL. 

j4[ 2 ^ N\ ^ 6 192 48iV 



+ 



For checking the convergence at NNLO, we carry out fits employing both M^/f'^ and M^/ 
as the expansion parameter. For the former fit we use (j30p to reexpand the expression for 
the form factor to the necessary order, and use the value of / obtained in [10]. Besides the 
pion decay constant the ChPT formula of the form factor to NNLO depends on 5 other LECs: 
II — ^2/2, ^4, ^6' ^yi' very difficult to find a stable fit in the 5-dimension parameter space. 

Therefore we fix li,l2 at the phenomenology values since they were calculated with small error 



bar from experimental data 



24( 1 ■ For which is only required in the formulation with f.,^, we 



apply the value obtained by an NLO fit of data in the range < 411 MeV by the PACS-CS 
collaboration [lo| ■ 

We find that stable fits are difficult to obtain unless we utilize data at more than a single pion 
mass in the fit procedure. The fit result obtained with data for the two pion masses = 296 
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4 G data at M =41 IMeV 
NNLO fit using f at M =41 IMeV 
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(a)NNLO SU(2) ChPT fit of tlie form factor utilizing /. {h){r^){fm^) calculated from the fit employing /. 

FIG. 9. Results of NNLO SU(2) ChPT analysis using the decay constant / in the SU(2) chiral hmit and 
combining data at M-^ — 296 MeV and 411 MeV. In panel (b) filled circles are values estimated from 
monopole fits, and open circles and lines are fit results. For the latter NLO and NNLO contributions are 
also plotted. 
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(a)NNLO SU(2) ChPT fit of the form factor utilizing (h){r'^){im^) calculated from the fit employing 

FIG. 10. Results of NNLO SU(2) ChPT analysis using the decay constant f^^ measured at each A/jr and 
combining data at M-^ = 296 MeV and 411 MeV. In panel (b) filled circles are estimations from monopole 
fits, and open circles are fit results. Upward triangles and downward triangles are contributions to NLO 
and NNLO contribution, respectively. 

and 411 MeV, which were used for the NLO analysis, is Usted in Table rvTIIl and shown in Figs. [9] 
and [TOl Combining the data at the two pion masses is acceptable since strange quark mass does 
not vary much, rus = 89.8(1.3) and 92.2(1.3) MeV at = 296 and 411 MeV, respectively 0. 
We observe that the results for E are consistent between the two fits within the error of 10-15% 
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and so are the predictions for the squared charge radius at the physical point, indicating that 
ChPT reasonably converges at NNLO up to ^ 400 MeV and ?s 0.01 GeV^. This is also 
seen by plotting the NLO and NNLO contributions separately as shown in Figs. [9] and [lOl The 
squared charge radius predicted at the physical point is close to the experimental value and is 
consistent within statistics errors of 10%. 

V. TOWARD THE PHYSICAL POINT A 64* LATTICE CALCULATION 

While ChPT to NNLO yields a reasonable result for the physical pion charge radius, the 
estimated error of 10% is quite large. We feel that for a convincing understanding of the 
physical pion charge radius one needs to explore the region of pion mass closer to the physical 
point than the value « 300 MeV analyzed so far. 

The PACS-CS gauge configurations has one more set corresponding to Mj^ ~ 156 MeV. We 
tried to calculate the form factor on this set, and found that the pion two- and three-point 
correlators exhibit very large fluctuations, to the extent that taking a meaningful statistical 
average is difficult. This trend becomes more pronounced as the twist carried by quarks becomes 
larger. Since LAIt^ ~ 2.3 at this pion mass for L = 32, we suspect that this phenomenon is 
caused by a small size of the lattice relative to the pion mass, and consequent increase of large 
fluctuations. 

A natural remedy to this difficulty is to employ larger lattices as one moves toward the 
physical pion mass. PACS-CS collaboration has been pushing a simulation on the physical point 
on a 64"^ lattice as a continuation of the work on a 32^ x 64 lattice. The hopping parameter of 
the run is adjusted to the best estimate of the physical point (k^,^. Kg ) = (0.137785,0.13665). 
We have used a subset of those configurations to calculate the pion form factor on a 64^ lattice. 
This requires much computer time, and hence we have only 4 configurations measured so far. 
We used the same setup as for the data set II, namely, (i) the incoming and outgoing pions in 
the three-point function carry momenta of the same magnitude \p'\ = \p\ but point in different 
directions, (ii) 4 random Z{2) Z(2) wall sources located at t = 0, 16,32,48 are employed, (iii) 
the twist technique is applied to the two quarks running from the source to the current and 
from the current to the sink, and (iv) four values are chosen for the twist angle 9 = (9, 0, 0) 
and its permutations such that the four-momentum transfer of the current takes the value 
Q2(GeV2) = 0.02,0.04,0.06,0.08. The fixed sink time tf is chosen to be 28. 
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FIG. 11. Triangle symbols represent effective mass obtained from measurement of 4 configurations of a size 
64^ with 4 source points located at t = 0, 16, 32, 48 and 4 random sources for each source location. Center 
lines exhibit PACS-CS estimation of pion mass calculated from larger statistics of 53 configurations. 

Since we use 4 random sources for each of the 4 locations of the source in time, our mea- 
surement on 4 configurations gives 64 two-point functions. The pion effective mass from our 
measurement together with the PACS-CS estimate of pion mass calculated from larger statistics 
of 53 configurations is shown in Fig. [TTl Although only 4 configurations have been used, one 
can already observe a plateau-like behavior for pion effective mass in this figure. The central 
value from the PACS-CS estimate corresponds to w 135 MeV. This is somewhat smaller 
than the charged pion mass, and significantly smaller than M-,^ = 156 MeV considered earlier 
on a 32^ x 64 lattice where we encountered problem of large fluctuations. 
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(a)Ratio i?'(r) defined in © at Q^(GeV^) = 0.02, 0.04. (b)Ratio R' {t) defined in (O at = 0.06 GeV^ 

FIG. 12. Ratio for extracting the form factor at several values measured on 4 configurations of a 64'* 
lattice. Fitting range is chosen to be t = 12 — 16. 
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FIG. 13. Momentum dependence of the form factor at ~ 135 McV with measurement taken on 4 
configuratfons of a 64'' lattice. Black line is a fit of two data points closest to = to the monopole 
ansatz (|2T|) . 

In Fig. [12] we plot the ratio R'{t) obtained from the 4 configurations. Making a constant 
fit over r = 12 — 16 yields the result for the form factor displayed in Fig. [131 Our data seems 
reasonable up to = 0.06 GeV^. Estimating the slope at the origin by a monopole fit to the 
two points closest to = 0, we obtain (r^) = 0.675(285)fm^. While the error is too large to 
seriously discuss consistency with experiment, it is certainly encouraging that the value is larger 
than those obtained at of about 300 MeV, and that the physical point simulation appears 
possible for the pion electromagnetic form factor on a 64^ lattice. 

VI. CONCLUSION 

We have presented a lattice calculation of the pion electromagnetic form factor in 2+1 dy- 
namical flavor QCD with the 0(a)-improved Wilson-clover quark action and Iwasaki gauge 
action. 

In order to obtain data with reasonable error at light quark masses close to the physical 
point, we have utilized some improved techniques besides traditional methods for the form 
factor calculation. We have shown that, choosing momenta of the incoming and outgoing pions 
to have the same magnitude but different directions, the ratio for extracting the pion form 
factor becomes statistically much better behaved. We have confirmed the validity of the twisted 
boundary condition and employed it to explore the form factor in the region of small four- 
momentum transfer. Application of the random Z{2) (8) Z{2) wall source has helped us to save 
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computing time considerably. 
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FIG. 14. (r^) (fm^) in comparison witli previous studies. The left-most filled circle represents our NNLO 
SU(2) ChPT prediction at the physical pion mass. 

Our data for the pion mean-square charge radius (r^) agree with recent data of other groups, 
and show that (r'^) increases toward the physical value as M^^ decreases. Nevertheless, on a 
32^ X 64 lattice, we could extract reasonable data only down to k, 296 MeV. 

ChPT analysis of the form factor utilizing NLO SU(3) or SU(2) formula lead to the squared 
charge radius which is smaller than experiment by 15-20%. Employing ChPT to NNLO improves 
the agreement. In fact our NNLO SU(2) fit using M^/f^ as the expansion parameter yields 
{r'^)pi^yg = 0.441(44)(13)(fm^) where the first error is statistical and the second error due to the 
error in the lattice spacing. 

We feel that a complete explanation of the behavior of the squared charge radius would 
require successful calculation of the form factor below ~ 300 MeV. Our experience points 
toward the necessity of enlarging the lattice size sufficiently. A sample calculation on a 64^ 
lattice indicates that the 64^ lattice with M^^L ~ 4 probably satisfies the requirement. We leave 
further exploration of the form factor calculation on such a lattice as future work. 
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